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Energetic particle irradiation of solids can cause surface ultra-smoothening ^, self- 
organized nanoscale pattern formation ^, or degradation of the structural integrity 
of nuclear reactor components ^. Periodic patterns including high-aspect ratio quan- 
tum dots with occasional long-range order ^ and characteristic spacing as small 
as 7 nm have stimulated interest in this method as a means of sub-lithographic 
nanofabrication ^. Despite intensive research there is little fundamental understand- 
ing of the mechanisms governing the selection of smooth or patterned surfaces, and 
precisely which physical effects cause observed transitions between different regimes 
1^ has remained a matter of speculation Here we report the first prediction of 
the mechanism governing the transition from corrugated surfaces to flatness, using 
only parameter-free molecular dynamics simulations of single-ion impact induced 
crater formation as input into a multi-scale analysis, and showing good agreement 
with experiment. Our results overturn the paradigm attributing these phenomena 
to the removal of target atoms via sputter erosion. Instead, the mechanism dominat- 
ing both stability and instability is shown to be the impact-induced redistribution of 
target atoms that are not sputtered away, with erosive effects being essentially irrel- 
evant. The predictions are relevant in the context of tungsten plasma-facing fusion 
reactor walls which, despite a sputter erosion rate that is essentially zero, develop, 
under some conditions, a mysterious nanoscale topography leading to surface degra- 
dation. Our results suggest that degradation processes originating in impact-induced 
target atom redistribution effects may be important, and hence that an extremely low 
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sputter erosion rate is an insufficient design criterion for morphologically stable solid 
surfaces under energetic particle irradiation. 

At irradiation energies between 10^ eV — 10^ eV, the irradiation process is domi- 
nated by the nuclear collision cascade caused by ion impact 03111. Displaced atoms that 
reach the surface with enough kinetic energy to leave are permanently sputtered away; 
all other displaced atoms come to rest within the solid or on the surface after phonon 
emission times of ~ 10~^^ seconds. These processes contribute prompt ero5/ve EMi and 
prompt redistributive \^^ \ ^^ \ ^^\ components of morphology evolution, respectively, and are 
collectively denoted P [x]. For most materials other than elemental metals, the damage 
resulting from these collisions quickly (~ 10^'^ seconds) leads to the amorphization of a 
thin layer of target material. Over much longer time scales (~ 100 seconds), mass trans- 
port by kinetic relaxation processes causes a gradual relaxational effect Ginzl|t] Hence, to 
the prompt term P [x] we add a phenomenological term for the gradual relaxation regime 
denoted G [x], assuming a mechanism of ion-enhanced viscous flow, which is expected 
to predominate in irradiated amorphous materials near room temperature . The prompt 
and gradual contributions to the rate of motion of the surface in the normal direction v-n 
superpose: 

v^ = P[^] + G[^]. (1) 

The prompt regime may be characterized using molecular dynamics (MD) simu- 
lations EHH or experimental methods Given data from many impact events, we may 
obtain the "crater function" A/i (x — x', 0) describing the average change in local surface 
height at a point x resulting from a single-ion impact at x', with incidence angle ^^j^ We 
then upscale the crater function into a continuum partial differential equation for the sur- 
face evolution using a multi-scale framework. The theoretical formalism for this process 
is described elsewhere here we provide a brief summary of the important points for 

^ Note that our focus on amorphous materials precludes potentially confounding effects due to singular 
crystallographic energetics and kinetics, such as the Villain instability in which surface diffusion is 
destabilizing and may be responsible for pattern formation on faceted single-crystal metal surfaces at certain 
temperatures'^. 

^ In principle, A/i also contains an explicit dependence on the initial surface shape /i(x) (as opposed 
to the implicit dependence that arises via the angle-dependence). This effect is described in but is very 
difficult to capture with molecular dynamics, so we leave its analysis for future work. 
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the linear case. Given the crater function A/i and the flux distribution / (x), we write the 
prompt contribution to surface evolution as a flux-weighted integral of the crater function 



A well-known observation in the field is that scale of the craters is much smaller than the 
scale of the resulting pattern and of the flux distribution. To exploit this fact, we employ a 
formal multiple-scale analysis, based on a small parameter e related to the ratio of impact 
scales to pattern scales. This formalism allows ready separation of the spatial dependence 
of the crater function (fast) from that of the surface shape (slow), and leads eventually to 
an upscaled description of the prompt regime of the form 



Here the Vs represent surface divergences, and the (increasing-order) tensors M^*) are 
simply the angle-dependent moments of the crater function A/i. This compact formula- 
tion is interesting for two reasons. First, the moments are readily obtainable directly via 
MD simulation, and they converge with far fewer trials than do descriptions of the entire 
crater function. Second, while atomistic methods have been used in the past to obtain the 
amplitude of a single term in a PDE obtained via phenomenological modeling I ^^ l ^^l t^, we 
believe this is the first derivation of an entire PDE from molecular dynamics results. 

Equation ([3]) describes the prompt regime P [x]. To fully capture the surface dynam- 
ics, we add to this a relaxation mechanism G [x] associated with ion-enhanced viscous 
flow Together, P [x] and G [x] completely determine surface morphology evolution 
via Equation ([T]). From this (nonlinear) equation, pattern-forming predictions are then ob- 
tained by examining stability of the the linearized equation as a function of the laboratory 
incidence angle 6. The derivation follows that in'^, and in an appropriate frame of refer- 
ence one finds that the magnitude of infinitesimal perturbations h away from a flat surface 
evolve, to leading order in e, according to the PDE 
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where the angle-dependent coefficients 

Sx{^) = Io^ [mW (0)cos(0)] 

. (5) 

Sy (0) = JoM^i) (0) cos (0) cot (0) 

are determined from the first moments obtained via MD, and the constant coefficient B 
is estimated from independent experiments. The structure of Equation indicates that 
linear stability is determined strictly by the signs of the calculated coefficients (Sx, Sy)' 
for values of 9 where either of these coefficients is negative, linearly unstable modes exist 
and we expect patterns, whereas for values of 9 where they are both positive we expect 
flat, stable surfaces. 

Existing uses of MD crater data for investigations of surface pattern-forming are en- 
tirely numerical in nature and could be viewed as a scheme for numerically integrating 
Equation |2j In contrast, our analytical upscaling of Equation |2] illustrates exactly which 
qualities of the crater - namely, its moments - play the dominant role in surface evolu- 
tion. Furthermore, this analytical form can be linearized, allowing predictions of stability 
boundaries, and changes to those boundaries as crater shape is varied. A crucial com- 
ponent of our approach is that the crater function Ah - and hence the moments M*^*) - 
contains the contributions of both erosion and mass redistribution. Whereas these effects 
have traditionally been treated separately via unrelated phenomenological models, view- 
ing the crater function as fundamental integrates erosion and redistribution into a unified 
description, allowing both processes to be treated identically and readily separated and 
compared. Indeed, this approach has permitted us to confirm for the first time conjectures 
E^mEHthat the stability of irradiated surfaces could be dominated by redistributive effects, 
with erosion - long assumed to be the source of roughening - being essentially irrelevant. 

We obtain angle-dependent moments from a series of MD simulations, in an envi- 
ronment consisting of an amorphous, 20x20x10 nm Si target consisting of 219,488 atoms 
created using the WootenAVinerAVeaire (WWW) method and then annealed with the 
EDIP potential This gives an optimized amorphous structure where most of the Si 
atoms have coordination number 4. The target was then bombarded with Ar at 100 eV and 
250 eV. During bombardment, the interaction between Si atoms was again described using 
the EDIP potential, which gives a good agreement between simulated and experimental 
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sputtering yields , whereas the Ar-Si interaction was a potential calculated for the Ar-Si 
dimer The kinetic energy was gradually removed during the simulations from the 1 nm 
borders of the substrate to prevent it re-entering the impact area via the periodic bound- 
ary conditions used in the simulations. The ambient temperature in the simulations was 
K. The simulation arrangements and their suitability for cluster and ion bombardment 
simulations are discussed in more detail in the supplement and in Refs. HHHIHll]. 

For each energy, two hundred impacts were simulated at each incidence angle in 5- 
degree increments between and 90 degrees, yielding moments as summarized in Figure 
1 . From the initial and final atomic positions, moments were obtained using the method 
described in"^. Briefly, by assuming that densities in the amorphous layer attain a steady 
state (i.e., that defect distributions immediately project to the target surface), we obtain 
erosive moments by assigning a height loss at each location proportional to the number 
of sputtered atoms originating from that location, while redistributive moments were ob- 
tained by assigning height losses at initial atomic positions and height gains at final atomic 
positions. For use within our analytical framework, we also fit the moments to Fourier se- 
ries constrained by symmetry conditions and by the observation that all moments tend to 
zero at 6^ = ±90°. For both energies, the redistributive first moments are much larger in 
magnitude than - and have the opposite sign of - the erosive moments. The implication 
is that redistributive effects completely dominate erosive effects, except possibly at the 
highest (grazing) angles where all moments tend to zero. 

To corroborate this finding, we calculate in Fig. 2 the coefficients in equation ([5]) for 
the 250 eV moments, and compare the pattern wavelengths they predict to experimental 
observations in the same environmental conditions ^'^ (clean linear experimental data at 100 
eV are not currently available). The agreement is generally good at intermediate angles, 
but several discrepancies between theory and experiment should be addressed. First, the 
small quantitative difference between predicted and observed bifurcation angles - which 
depends only on the shape of Sx (d) - could readily arise from the approximate nature 
of the classical potential, on which our simulations are based (unlike the bifurcation an- 
gle, precise wavelength values depend on the value of B, which could only be estimated). 
Second, the measured moments do not predict a transition to perpendicular modes at the 
highest angles; this could be due to our neglect of explicit curvature-dependence in the 
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crater function, but additional physical effects such as shadowing and surface channeling 
- not addressed here - are known to be important at grazing angles. Third, we find no pre- 
diction via MD of the experimentally-observed perpendicular-mode ripples at low angles-; 
indeed, our redistribution-dominated continuum PDE is maximally stable at low angles, 
and equations of the general form|4]are anyway unable to generate the constant- wavelength 
or "Type I" bifurcation that is observed 'i^. These observations, together with the experi- 
mental observation that low-angle ripples develop over much longer timescales than their 
high-angle counterparts '23, suggest that low-angle perpendicular- mode ripples are not due 
to crater-function effects at all. It has already been observed ^° that the low-angle Type I 
bifurcation is consistent with any of several non-local physical mechanisms such as grad- 
ual stress buildup and relaxation or non-local damping None of these effects would 
be captured in the (prompt, local) crater function, and this observation motivates future 
studies incorporating such physics. 

Despite the limitations of our approach, when one considers the lack of any free pa- 
rameters in the theory, the agreement for the diverging-wavelength or "Type 11" bifurcation 
^ near 50 degrees is remarkably good. The agreement remains good even when the ero- 
sive coefficients are omitted, and the similar shapes of the redistributive moments at 100 
and 250 eV is consistent with the reported ^ energy-insensitivity of the stability bound- 
ary. The most striking aspect of this result is its logical conclusion that erosion effects are 
essentially irrelevant for determining the patterns: according to Fig. 2 the contributions 
of redistributive effects to the S coefficients, which determine stability and patterns, are 
about an order of magnitude greater and opposite in sign. This conclusion overturns the 
erosion-based paradigm that has dominated the field for two decades and we suggest 
its replacement with a redistribution-based paradigm. An important direction for future 
research is identifying the energy range over which this conclusion holds. Although it is 
conceivable that erosive effects might become non-negligible at higher ion energies, or 
for dense metals where heat spikes enhance sputtering yields, it remains to be determined 
whether erosion is actually important for stability or pattern formation in any physical 
experiment to date. Preliminary analysis of simulations at 1 keV are consistent with our 
conclusions from the results at 100 eV and 250 eV, and experimental results at 1 keV lead 
to the same conclusion!^. 
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Non-erosive ion impact-induced atom redistribution at surfaces has heretofore not 
been considered in the design of plasma-facing fusion reactor wall materials, where low 
sputter yield has been an important design criterion in the selection of tungsten for sta- 
ble surfaces that must be exposed to large plasma particle fluxes for extended periods. 
Because the average helium ion energy is only ~ 60 eV and the threshold energy for 
sputter removal of tungsten is ~ 100 eV, this material has been considered impervious 
to the effects of erosion. Within the erosion-based paradigm of pattern formation, the 
nanoscopic surface morphology that evolves on tungsten surfaces under some such condi- 
tions has therefore appeared mysterious ^. However, because the negligibility of erosive 
effects does not prevent redistributive effects from causing pattern-forming instabilities, 
as we have shown here through a crater-function analysis, atom redistributive effects may 
be important contributors to the origin of these mysterious morphologies, and we present 
in the Supplement a re-analysis of data gathered in 1^ that supports this idea. If this con- 
jecture turns out to be correct, then ultimately crater function engineering considerations 
may provide a more refined materials design criterion than simply a low sputter yield. 
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Figure 1: Fitted angle-dependent moments of the crater function Ah (x, y) as determined 
from molecular dynamics, for both 100 eV and 250 eV. (a) Zeroeth erosive moment M*^°^ 
(sputter yield times atomic volume). (b,c) First erosive (b) and redistributive (c) moments 
M^^\ Each of the latter contain components in both the x (downbeam) and y (crossbeam) 
directions, with the latter expected to be zero from symmetry arguments At both ener- 
gies, redistribution dominates erosion. 
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Figure 2: Comparison between predicted and measured wavelength for 250 eV Ar— )■ Si. 
Dimensional coefficients were calculated with a flux of / = 3.5 x 10^^ ions/ (cm^s) for 
comparison with experiment'^. (a,b) Coefficients of |^ and |^ in the linearized evolution 
equation|4]using the experimental flux. These coefficients are dominated by redistributive 
effects, (c) Comparison of predicted ripple wavelengths with average experimentally- 
observed wavelengths. Circles / squares indicate experimental patterns with wavevector 
parallel / perpendicular to the ion beam, and the vertical dashed black lines indicate experi- 
mental phase boundaries. On top of this, the solid line indicates our predicted wavelengths 
(which are all parallel-mode), and the dashed blue line indicates the wavelengths predicted 
if erosion were neglected entirely. 
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A Simulation Methods 



In the simulations, the actual impact induced structural deformations can be detected only 
in a well relaxed silicon structure. A stress existing within the material as a result of incom- 
plete relaxation of the structure before the simulation can induce displacements of silicon 
atoms upon impact. For example, our test simulations showed that a structure created by 
annealing sihcon in MD is not dense enough to model the real amorphous silicon. The sur- 
face of such a structure collapses upon impact. It is possible to relax the internal stresses by 
bombarding the surface with silicon atoms before the actual simulation. However, the re- 
laxation will be not uniform. Therefore, we have used the WootenAVinerAVeaire (WWW) 
methodic. 

The WWW method is a computer algorithm that generates realistic random network 
models of amorphous silicon. In this method, the structure is described with positions of 
N atoms and a list of bonds between the atoms. After a random switch of two bonds, 
the structure is relaxed using an interatomic potential HIMI. in connection to the struc- 
tural relaxation, the Berendsen pressure control algorithm is used to relax the diagonal 
components of the pressure tensor to zero The algorithm is computationally demand- 
ing and therefore it is possible to fully optimize only relatively small structures of a few 
thousand atoms. Therefore, the optimized block must be copied and the copies joined 
to create an optimized silicon structure large enough for impact simulations. Our tests 
showed that the best approach is to partially optimize a rather large (10x10x10 nm) build- 
ing block instead of building the structure of very small fully optimized blocks. The latter 
approach can induce artificial internal shear stresses in spite of the optimization (note that 
pressure control at periodic boundaries an an orthogonal cell does not necessarily remove 
the nondiagonal, shear components of the stress tensor). The optimization of a 10x10x10 
nm amorphous silicon structure was achieved using a parallelized implementation of the 
WWW algorithm. 

After the WWW optimization phase, the structure and the surface were relaxed in 
MD using the EDIP potential, before it was used in the actual impact simulations. The 
structure used in the simulations was built of four identical optimized blocks. The den- 
sity was 2.5g/cm^, which indicates that the structure is dense and not likely to collapse 
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upon impact. Two-thirds of the silicon atoms had four neighbours, the others had five 
neighbours, which is also a sign of a dense structure. In a perfect amorphous network, 
all silicon atoms should have four neighbors. However, the test with structures built of 
smaller blocks that were better optimized than the 10x10x10 nm structure showed that 
the behavior of the moments as a function of the impact angle are very similar as those 
reported in the main article. 

To simulate ion irradiation of bulk silicon samples, we have used periodic boundary 
conditions with slab boundary conditions (free boundary in the incoming ion direction, 
and periodic boundaries in the other two lateral directions). The bottom 1 nm layer of 
atoms in the simulation cell were held fixed, and temperature scaling was also applied to 
the atoms in a 1 nm thick layer above it. A border region of 1 nm thickness in the lateral 
directions was cooled during the simulations. The role of the cooling zones is to prevent 
shock waves and phonons to re-enter the impact area via the periodic boundaries. The size 
of the simulation box (20x20x10 nm) was chosen to be large enough to contain not only 
the area containing the crater but also the area of small deformations which reach 5-7 nm 
from the impact point. 

In addition to the good amorphous silicon structure and cooling of the boundaries, 
the quality of the repulsive Ar-Si potential affects the outcome of impact simulations. 
The kinetic energy of the impacting Ar atom is deposited first in relatively strong Ar-Si 
collisions where the atoms come close each other. The repulsive Ar-Si potential used in 
this study is calculated using density-functional-theory calculations utilizing a numerical 
basis sets. This method is shown to give a more accurate repulsive potential than the 
standard universal ZBL potential The same method was used to create Si-Si short-range 
repulsive potential which was smoothly joined to the EDIP potential. This ensures that 
the possible collisions between high-energy primary knock-on silicon atoms are correctly 
modelled. 

The same initial structure was used for all simulations. However, to reduce as much 
as possible the effect of the initial surface structure on the measured crater statistics, the 
following steps were taken. First, the impact point was varied randomly over the entire 
surface, with (periodic) cooling regions dynamically re-assigned for each simulation so as 
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to be maximally distant from the impact point. Second, the azimuthal angle was varied 
randomly, and only the atoms within 7.5 nm of the impact point were included in the 
moment calculations (i.e., cooling zones and corners of the square box were ignored). 
Last, the global average - or "background" - displacement (over all simulations at a given 
impact angle) was subtracted from each individual displacement before analysis. Because 
we varied the impact point and azimuthal angle randomly, this global average ought to 
be zero for a perfect target over enough trials, and subtracting it helps to remove target- 
specific effects from our measurements. With this setup, 200 simulations at each impact 
angle were performed, which was sufficient to isolate the background displacement, and 
to cancel the effect of thermal vibrations which were present in the structure after impact. 

A final challenge in analyzing the data arose due to a combination of the amorphous 
nature of the target with the periodic boundary conditions. On an ideal, very large MD 
target, the effect of the ion impact would only be felt within a finite distance from the 
impact point. However, to allow the gathering of data within reasonable time, a limited 
box size must be used, and the periodic box described above - with cooling zones - appears 
to be the most physically plausible way of accomplishing this. However, a periodic box 
always means that one is truly simulating an infinite number of simultaneous side-by-side 
impacts, and with a small enough domain, these impacts can generate enough co-ordinated 
momentum transfer to shear the entire target, especially for an amorphous target. Indeed, 
within our target, the average downbeam displacement of atoms was consistently a linear 
function of distance from the target's rigid floor, except near the surface, where larger 
displacements were concentrated. 

To combat the problem of "global shear," we measured the average downbeam dis- 
placement within different annular slices of target parallel to the surface, using inner/outer 
radii of 2nm/9nm (i.e., away from both the impact and the cooling boundaries). We then 
fit the bottom half of the resulting depth-dependent profile to a line using least squares, 
and finally subtracted the extrapolation to the surface of this fit from the overall displace- 
ment field, as illustrated in Figure[3] The results localize the displacements to within a few 
nanometers of the surface, which is consistent with measured amorphous layer thicknesses 
of 3 nm for Si irradiated by Ar at 250 eV'^. In the future, we will explore the response of 
a larger, hybrid target consisting of a 3 nm layer of amorphous Si atop a crystalline (but 
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not rigid) base. However, we believe our existing measurements are within the accuracy 
level of the other estimates in the paper. 



B Moment Capture and Fitting 

Here we describe in more detail how we obtain moments, how we fit them, and how we 
obtain the final linearized evolution equation (4). For each simulated ion impact with our 
initially flat MD target, we define a co-ordinate system (x, y, z) centered at the impact 
point with z pointing normal to the surface, x pointed in the direction of the projected ion 
path, and y perpendicular to both x and z so as to complete a right-handed co-ordinate 
system. Hence, the ion is always arriving from the negative x-direction. Hence, the crater 
function Ah (x, y; 6) describes the height change associated with an impact at the origin, 
of an ion with indicence angle of 6 from the vertical, coming from the negative x-direction. 

After impact, we extract the moments 



using the method described in Here, it is important to note that the first moments 
M^^) contain both erosive and redistributive components (because the zeroeth moment 
M^'') describes mass loss, and because redistribution is mass-conserving, M*^°) has no 
redistributive component). 

Simulation sets were performed at 5-degree increments, and the average of the re- 
sulting moments were fit to Fourier functions of the form 




Ah (x, y) dx dy 




(6) 




Ah (x, y) y dx dy 
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3 

M^_odd = '^o.n sin {2n9) 

V (7) 

M^_even = &n COS {{2n - 1) 9) 

n=\ 

These fittings reflect the observation that aU moments tend to zero at 6' = 90°, and also 
their symmetries about 6' = (because a positive theta indicates an ion beam coming from 
the negative x-direction, moments that are odd/even in x should also be odd/even in 6'). 
As seen in the main text, this method does not produce perfect fits, but the use of simple 
Fourier modes eliminates potential model-bias, while the restriction to a small number of 
terms excludes inter-angle noise from the fitted curve. 

For our data, the moment fits were given by: 

M(°) ^ -12.0 cos (Q) + 13.2 cos (3^) - 3.46 cos (5^) [lO^^ nmVion] 
^iiLive ~ -10.4 sin (2^) + 5.81 sin (4^) - 0.516 sin (6^) [lO~^ nmVion] . (8) 
^Sdist. ~ 291 sin (2^) - 39.6 sin (40) - 2. 16 sin (60) [lO"^ nmVion] 

The moments m'^^ are zero to within sampling error, as expected from symmetry consid- 
erations. 

C Analysis: from Moments to Coefficients 

For the general reduction of moments to (nonlinear) PDF terms P [x], we refer to the 
framework derived in However, in the linear case discussed here, it is sufficient to 
consider the linearization of the first-order term obtained by combining Fquations (1) and 
(3) from the main text: 

vl (x) ^ £Vs ■ (/M(i)) (9) 

(in the linearization, the zeroeth- and second-moment terms do not contribute to stability). 
Now, V5 indicates a surface divergence, and indeed this calculation is most naturally 
performed in a local co-ordinate system associated with the surface normal and projected 
beam direction. In particular, both the flux / (0) and the moments M*^*) (0) depend on 
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the local incidence angle 0, while the vector M'^^^ (0) is observed to always point in the 
direction ep of the projected ion beam. 

Following surface velocities at an arbitrary point x will be calculated in a local 
co-ordinate system (U, V, W) centered at x, where e^y = n corresponds to the surface 
normal, eu = ep corresponds to the downbeam direction associated with the projected ion 
beam, and ey = e^y x e^/. In this system the surface is described locally by the equation 
W = H {U,V), the ion flux is / = Jq cos (0), and the first moment is M^^^ = f (0) ep, 
where / (0) = M^^ (6) as measured from MD. Now, as described in"^, it is sufficient for 
the purposes of calculating one surface divergence to approximate H {U,V) via 

H^^ [HuuU^ + 2HuvUV + HyvV^) . (10) 

where Huu, Huv, and Hyy describe the surface curvature at x. All other variable quan- 
tities can then be approximated in the vicinity o/x to first order in(f/, V) via: 

dH dH 

n 



dU' dV 

I ,dH dH\ 

ep«^l,-cot(^„)-.-^ (11) 

dH 

cos (0) ^ COS (0o) + — sin (0o) 

When we now take the surface divergence = {du, dy) and evaluate at {U, V) = 
(i.e., at x), we obtain in the local co-ordinate system 



where 



< (x) ^ eVs ■ (/M(i)) = elo [Su (0) Huu + Sy (0) Hyy] (12) 



^.(0)-|[/(0)cos(0)] _ ^^^^ 

Sy (0) = /(0)cOs(0)cot (0) 



While linear in the local co-ordinates. Equation ( [T2| ) is in general nonlinear in the lab fram. 
However, for stability studies we need only the linearization of ([12]) in the lab frame, which 
in dimensional form is simply 

dh 



dt 



prompt 



^o{Sx{e)^,+Sy{e)^] (14) 
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because, to linear order, 



Sx{e) 
SY{e) 



Hi 



uu 



vv 



(15) 



Sv (0) 



To expression ( [T4| ) for the prompt regime we add the linearization of the gradual regime 
associated with ion enhanced viscous flow which is a lubrication approximation with 
the form 

§ = -B^'h. (16) 

^'^ gradual 

Adding the prompt and gradual regimes, we obtain the evolution equation (4) in the main 
text. 



For completeness, we conclude with the functional forms of the coefficients {Sx, Sy) 



associatd with our fittings, which are obtained from equations ( 13 1: 



Qeros. 

ceros. 

oredist. 
"^X 

Qredist. 



-52.1 COS {6) - 69.2 COS {W) + 132 cos (5^) - 18.1 cos (7^) 
-50.5 cos {0) + 24.7 cos (3^) + 21.3 cos (5^) - 2.58 cos (7^) 
1450 cos {6) + 3760 cos (3^) - 1040 cos (50) - 75.6 cos (70) ' 
3520 cos (0) + 815 cos (30) - 231 cos (50) - 10.8 cos (70) 



(17) 



D Estimation of Viscous Flow coefficient 

The materials parameter B appearing in Equation (4) of the main text is defined!^ as 

= (18) 

where 7 is the surface free energy, d is the thickness of the thin amorphous layer that 
is being stimulated by the ion irradiation, and 77 is the layer's viscosity. We assume the 
surface free energy of amorphous silicon under ion irradiation to be equal to its value in 
the absence of irradiation; the value of 7 = 1.36 J/m^ measured via molecular dynamics 
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simulations by Vauth and Mayr ^ happens to be numerically equal to that measured ex- 
perimentally for single-crystal Si(OOl) For the amorphous layer thickness, we directly 
measured d ~ 3.0 nm via cross-sectional transmission electron microscopy on samples 
irradiated at normal incidence and 30 degrees from normal. Finally, we estimate the vis- 
cosity of the top amorphous Si layer during irradiation to be ^ 6.2 x 10^ Pa sec, as 
shown below. 



The reciprocal of viscosity is the fluidity 0, which is generally understood to scale 
with the flux / , and can be expressed in the form 

4> = Hx iVopAPs (19) 

where H is the radiation-induced fluidity, and A'^dpaps is the average number of displace- 
ments per atom per second. Using molecular dynamics simulations, Vauth and Mayr ^ 
report H = 1.04 x 10^^ (Padpa)^^ at an energy E = 1 keV and temperature T = 300K 
- we use this value for our comparison with experiment, with the caveats discussed below. 
The average number of displacements per atom per second is given by 

-^DPAPS = -^recoils 5 (20) 

a 

where O, = .02nm^/atom is the atomic volume of silicon, / = 3.5 x 10^^ ions/ (cm^s) 
is the experimental flux in the plane perpendicular to the ion beam, d ^ 3nm is the 
amorphous layer thickness, and iVrccoiis is the number of recoils generated per ion impact. 
To estimate A^iccoiisj we use the Kinchin-Pease model for the gross number of Frenkel 
pairs per incident ion, obtaining A^recoUs = O.SE/Ed = 6.7, where E = 250 eV is the ion 
beam energy and Ed = 15 eV is the displacement threshold energy of Si^. Taking all of 
these quantities, we obtain a value of the viscosity of 

r/=— ^ = 6.2 X lO^Pasec. (21) 

H X Adpaps 



As discussed above, the value for rj listed in ( [21] ) is associated with Vauth and Mayr's 
value of H = 1.04x 10^^ (Pa dpa)^^ at an energy E = 1 keV and temperature T = 300K. 
In contrast, our experiments were carried out at E = 250 eV, and the irradiated sample is 
observed to reach temperatures of approximately 450 K. Hence, there is some uncertainty 
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in our value of rj, which translates to uncertainty regarding the vertical position of the 
theoretical curve in Figure 2 of the main text. For the temperature difference, Vauth and 
Mayr observe H to increase weakly with increasing substrate temperature, which would 
shift the curve upward; however, the shift would likely be less than a factor of two. As 
for the energy difference, in a study of CuTi, another amorphous material, Mayr et al ^ 
found H to either increase or decrease with recoil energy depending on the details of the 
simulations; hence the theoretical curve in Fig. 2 of the main text could shift either upward 
or downward for MD simulations at 250 eV, with a potential magnitude of perhaps a factor 
of two. 

E Preliminary Supportive Data 

We believe the phenomena reported in this paper are applicable to a wide variety of sys- 
tems. However, the main data sets leading to the discovery of this result both describe 
low-energy irradiation of amorphous silicon. Therefore, we provide here some prelim- 
inary results associated with higher-energy argon irradiation of silicon, and also for the 
case of helium-irradiated tungsten which was referred to in our conclusion. 

For silicon irradiated by argon at 1 keV, a much larger target must be used, with 
dimensions of 40 x 40 x 10 nm, which makes simulations expensive. However, we have 
performed 30 simulations at 5-degree increments between 30 and 65 degrees, which spans 
the experimental transition, and the results so far are consistent with those at 100 and 250 
eV: the redistributive contributions to the first moment dominate the erosive contributions. 

For tungsten irradiated by helium at 100 eV, average displacements are so small that 
we were not able to isolate a clear signal against background noise after 200 simulations. 
However, using existing data from earlier workf^, we were able to calculate the cumula- 
tive displacement field at a single angle after 4,000 simulations. In all of these simulations, 
not a single tungsten atom was sputtered, yet a downbeam bias in the displacement field is 
clearly observed in Figure 3. Hence a crater- function approach may enable better engineer- 
ing of surface stability under conditions where the sputter yield is zero but impact-induced 
target atom displacements still occur. 



18 



2 



drift fit and correction 




""-0.1 0.0 0.1 0.2 0.3 0.4 0.5 

drift [10"'' nm] 



Figure 3: Illustration of the removal of shear at 60 degrees for irradiation at 250 eV. Green 
dots are the original measurements, and the blue line represents a linear fit to the bottom 
half of the dots, extrapolated to the surface. The red line is the result of subtracting this 
extrapolation from the original data, which localizes the displacements to the vicinity of 
the surface. All data are associated with a 2nm/9nm annulus that masks the impact zone 
and the cooling boundary zone. 



19 



W displacements by 4000 He ions at 100 eV ^ W 




'/y. 



Figure 4: Cumulative displacement field for 100 eV He — )• W after 4,000 impacts. Im- 
pinging helium ion comes from upper right, at an angle of 25 degrees from normal. The 
sticks combine the initial and final position of atoms that were initially in a half a unit cell 
thick cross section through the simulation cell. The displacements are analyzed for data 
initially simulated in Ref. The blue ends of the sticks indicate the initial and the red 
ends the final positions of the atoms. The surface is located at the top. No W atoms were 
sputtered, but a clear downbeam bias in the displacement field is seen. 
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